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Abstract 

We present a stochastic phase-field model for multicomponent lipid bilayers that explicitly ac- 
counts for the quasi-two-dimensional hydrodynamic environment unique to a thin fluid membrane 
immersed in aqueous solution. Dynamics over a wide range of length scales (from nanometers to 
microns) for durations up to seconds and longer are readily accessed and provide a direct compari- 
son to fluorescence microscopy measurements in ternary lipid/cholesterol mixtures. Simulations of 
phase separation kinetics agree with experiment and elucidate the importance of hydrodynamics 
in the coarsening process. 

PACS numbers: 87.16.A-, 87.16.dj, 83.10.Mj, 87.16.D-, 87.15.Zg 



Multicomponent lipid bilayer membranes are of universal biological importance and are 
increasingly viewed as fundamentally interesting soft matter systems [TH3]. Ternary mixtures 
of saturated and unsaturated lipids and cholesterol have become a standard experimental 
model to study the dynamics of inhomogeneous membranes under controlled laboratory 
conditions [HE]. These dynamics include: diffusion of lipid domains |6], "flickering" fluc- 
tuations of domain boundaries ^ and phase separation kinetics [8] . Analytical calculations 
for domain diffusion and flickering in certain regimes have helped to illuminate the rel- 
evant physics [H [10], which crucially depends on the "quasi-two-dimensional" (quasi-2D) 
[TT] hydrodynamic environment of a viscous fluid membrane suspended within bulk solvent 



Particle based simulations have been used to study general features related to phase 
separation dynamics in bilayer systems [HHTG], but it remains difficult to directly com- 
pare simulation to experiments that commonly involve vesicles tens of microns in diameter 
and observation times that run up to minutes. At a coarser level, phase-field representa- 
tions of inhomogeneous membranes have been introduced [TTHT^ . but these studies have 
neglected hydrodynamic effects and/or thermal fluctuations, again making detailed com- 
parison to experimental dynamics impossible. This letter presents a simulation scheme that 
combines quasi-2D hydrodynamics and thermal fluctuations with the capability to probe ex- 
perimentally relevant length and time scales. Compelling agreement with both theory and 
experiment is obtained, suggesting this methodology as a powerful tool to study membrane 
dynamics in the context of physics, biophysics and biology. 

Our focus is on experimentally observable dynamics, not first-principles prediction of 
detailed thermodynamic properties. Accordingly, we adopt a standard Landau- Ginzburg 
free energy functional for binary mixtures [201 12T] 



two-phase coexistence in ternary lipid/cholesterol systems. However, to facilitate comparison 
with experiment, we assume tight stoichiometric complexation between cholesterol and the 
saturated lipid[22] and identify 0(r) = unsaturated — x(i")compiex as the local difference in 
mole fraction (x) between lipid species. Though limited to mixtures with a 1:1 stoichometry 
between saturated lipids and cholesterol, this picture has the advantage of simplicity: only 



[T2l[T3]. 




(1) 



where r, u, 7 > 0. Eq. [T]may be interpreted as a phenomenological description of the observed 
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a single composition field is required and the three parameters appearing in F are readily 
related to measurable experimental properties [20]: the line tension between phases, a = 
^^^1/2^3/2^ interface width ^ = |a/7/t, and equilibrium phase compositions 0o = ±^/r/u. 

We require the dynamics of to conserve lipid concentrations, hydrodynamically couple 
points on the membrane and have thermal fluctuations. Model H dynamics represents the 
generic long- wavelength low-frequency picture to incorporate these requirements |2T]. In 
the overdamped "creeping-fiow" limit for experimental conditions (low Reynolds number, 
Re < 10~^), model H reduces to 



5F 

(5^ + V ■ V)0(r, t) = MV'—— + ^^(r, t) (2) 

d(p{v,t) 



^V;.0(r',t) + O(r',t) 



v,iv,t) = j d?r'Tij{Y-v' 

Though typically applied to pure three dimensional (3D) or two dimensional (2D) geometries, 
we use Eq. [2] for the quasi-2D fluid membrane geometry first introduced by Saffman and 
Delbriick |12], which considers the membrane to be a thin, flat fluid surface with surface 
viscosity rjm surrounded by a bulk fluid with viscosity rjf. In this case v = {vx.Vy) is 
the in-plane membrane velocity field, M is a transport coefficient related to the collective 
diffusion coefficient for lipids within the bilayer {D^) via M = D^/2r [HI |2T] and Tij[r) is 
the Green's function for in-plane velocity response of the membrane to a point force [TH |T3] . 
In a conventional 3D fluid geometry, Tjj(r) would be the Oseen tensor [25]; its form for 
the quasi-2D membrane geometry includes the effect of flow both within the membrane and 
in the bulk solvent. There is no simple closed form expression for Tjj(r), but its Fourier 
transform is [TTl [13] 

(q) = , ,j , (3) 

where Lsd = ^ is the Saffman-Delbriick length scale. The Gaussian white thermal noise 
terms, 9{r,t) and (j{r',t), are distributed with variances set by the fluctuation-dissipation 
theorem [26] . In Eq. [2] and henceforth the Einstein summation convention is assumed. 

In typical model membrane systems, 1]^ ranges from (0.1 — 10) x 10^^ surface poise 
(poise-cm, or grams/s) [SIEZIES], corresponding to Saffman-Delbriick lengths ~ 0.1 — 10 
microns. Eq. [3]exhibits the characteristic scaling (1/r in real space) associated with the usual 
3D Oseen tensor for q <^ L~]^, and reduces exactly to the 2D analog to the Oseen tensor for 
q ^ [llj. However, fluorescence microscopy experiments probe wavelengths comparable 
to and it is critical to use the full expression in comparison to these experiments. 
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Koga and Kawasaki [23] suggested the use of fast Fourier transforms as an efficient way 
to numerically evolve zero-temperature overdamped model H dynamics. We have extended 
their approach to include stochastic thermal forces and the quasi-2D hydrodynamics dis- 
cussed above, evolving Eq. |2]in Fourier space. Our evolution uses a Stratonovich scheme [30] 
with semi-implicit terms similar to those used for the deterministic Cahn-Hilliard equation 
[31]. This simulation methodology may be viewed as an extension of traditional "Brown- 
ian dynamics with hydrodynamic interactions" [32] to composition dynamics within a flat 
quasi-2D membrane environment. Details of the numerical scheme are in the appendix. 

The translational diffusion of circular lipid domains presents an ideal test to assess the 
validity of the simulation method. Theoretical results based upon the quasi-2D hydrody- 
namics underlying our approach have been derived [IDl US] and confirmed experimentally 
[6l[28l|33]. These results predict Da = {kBT/4:7irim)J-'{a/ Lgd) for the diffusion coefficient of 
a domain of radius a. The function J-" represents the solution to an integral equation [TO] , 
but is well approximated by a closed- form empirical fit described in [2SJ- Choosing initial 
conditions to refiect a single circular domain and thermodynamic parameters that guaran- 
tee the domain remains nearly circular (high a, low T), we track the domain's position over 
time and infer Da via its mean square displacement. The results collapse onto J^{a/ Lgd) (all 
of a, rim and rjf were independently varied) over a wide range of a/Lgd ratios, completely 
spanning the crossover between Saffman-Delbriick diffusion Da ~ ln(Lsd/a) in the limit of 
small domains and Stokes-Einstein-like diffusion Da ~ 1/a in the limit of large domains 

(Fig.g. 

In experimental membrane systems, line tensions are seldom high enough to fully sup- 
press shape fluctuations of lipid domains. Though these fluctuations have minimal impact on 
translational diffusion, the fluctuations themselves may be analyzed [7] and provide a further 
test of our simulation methods. Evolving an initially circular domain by Eqs. [2| we found, 
in qualitative agreement with [34J, that domains with small line tensions ~ 0.1 — 0.2 pN 
are not thermodynamically stable at temperatures ~ 20°C and area fractions ~ 0.03 — 0.1; 
the domain radii shrink over time in favor of a more homogeneous distribution through the 
simulation box (Fig. |2]). (However, this homogeneous phase does not appear to be composed 
of "an ensemble of small domains" [31]; the composition is rapidly fluctuating everywhere, 
without any clearly deflned domain boundaries.) This instability depends not only on line 
tension, but also on the lipid composition of the membrane; for 1:1:1 mixtures (area frac- 
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FIG. 1: (Color online). Translational diffusion coefficients obtained by tracking simulated domain 
motion, compared with Saffman-Delbriick-Hughes-Pailthorpe- White theory (using the interpola- 
tion of [28]). Dq = ksT /Aurira- Error bars are on the order of marker size. Line tensions and 
temperatures were chosen to ensure domains remained circular over the course of the simulations; 
box sizes range from 10 /xm to 40 /im (N = 256 to 1024), chosen to eliminate effects from periodic 
boundary conditions. 

tion 50%) with similar physical properties, micron-scale phase separation is observed both 
numerically and experimentally (Fig. |3]). 

At higher line tensions, domains under similar thermodynamic conditions are stable and 
we have simulated the dynamic fluctuations of these systems (Fig. |2| We use the image 
analysis techniques of [7j, tracing the boundary of the domain, and expanding it into quasi- 
circular modes, r{9,t) = Ro{l + Uo(t) + ^J2n=/=o'^ri(t)e^"'^). For small deviations m„, the 
domain shape is expected to behave in accord with an effective Hamiltonian H = aL ^ 
^^y^Z]n>o(^^ — 1)|m„P {L is the domain perimeter); thus (I-Unp) = 2kBT / a-nRoin^ — 1) is 
expected via equipartition. Experimentally, this relation is used to determine line tensions 
[7]. We observe that our numerical experiments yield the equipartition result with the 
expected a = ^^^^/^^^/^ (Fig. [2] inset). The linear dynamics of these modes depend on 
the line tension as well as the viscosities of both the membrane and the surrounding fluid; 
(M„(t)'u_„(0)) = (|M„p)e^*/'^" . Stone and McConnell found these relaxation times neglecting 
the membrane viscosity (valid for length scales R/n ^ Lgd) and Mann et al. calculated 
them neglecting the bulk fluid (valid for R/n ^ Lgd) [SI ES]- In our simulations, we see 
exponential relaxation over all modes u„(t) with x„ plotted in Fig. |2] We note that for 
large n modes, the Mann theory is appropriate, as expected, but there are deviations from 
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the Stone-McConnell result even at low n since Lgd = 0.5/im is comparable to the size of 
the domain. The simulations are in excellent agreement with a recent generalization to the 
Stone-McConnell theory that includes the effects of both rjm and rjf. For reference, we 
include the three expressions here: 



stone 



general 



VmRp 1 

cr n'^{n'^ — 1) 



dx 



na 



a;2(x + Ro/Lsd) 



where J„(a^) is a Bessel function of the 1st kind. 



(4) 
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FIG. 2: (Color online). TOP: Simulated evolution of a flickering domain with initial radius R = 
2.0/im, = 1 X 10^^ surface poise, a = 0.8,0.1 pN and ^ = 40 nm, 0o = 0.4. The collective 
diffusion coefficient = = 3.2 x 10"^ cm^/s. = 512, and the time step At is 5 /is. 

T = 21°C, rjf = 0.01 Poise. These are typical values for ternary domains near 20°C, with a 
viscosity in the middle of the commonly accepted range (see text). The cr = 0.1 pN domain 
shrinks, whereas the a = 0.8 pN domain is stable. The box size for the simulation is 20 ;um x 
20 jjLvn (Scale bar: 2 /im. Only the region immediately surrounding the domain is displayed.). 
BOTTOM: Relaxation times for mode n for a = 0.8 pN (see text). The theoretical predictions are 
summarized in Eq. [4} Inset: thermal variance in amplitude for mode n. The theoretical prediction 
(black) is the equipartition result described in the text. There is no fitting involved for either the 
main figure or inset. Theoretical predictions follow immediately from the parameters input to the 
model and demonstrate the robustness of the simulation scheme. 



The modeling described herein displays its true potential when applied to problems that 
are not easily explained with analytical theory. Using identical methodology to the diffusion 
and fluctuation studies, but employing a homogeneous initial condition, allows the study 
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of phase separation kinetics in ternary lipid/cholesterol model systems. Our simulations 
are motivated by the experiments of Veatch and Keller [H] on roughly 1:1:1 mixtures of 
DOPC/DPPC/Chol. Taking care to choose parameters consistent with the experimental 
system, we find strong qualitative similarity between simulation and experiment (Fig. |3]). 
However, we stress that this comparison has limitations. The experimental temperature 
quench taken in [S] was not precisely controlled or recorded; our homogeneous starting 
point, followed by constant T dynamics represents a numerically convenient choice, adopted 
in the absence of clear experimental guidance. A further limitation is that we have assumed 
both phases share the same viscosity. Lipid diffusion coefficients in ordered and disordered 
phases differ roughly by a factor of ten [39] , suggesting a comparable difference in viscosities 
between the two phases. The adoption of a single viscosity for the two phases is an approx- 
imation; nevertheless, certain composition dynamics in ternary model systems appear to be 
adequately described using single-viscosity theories employing (either explicitly or implic- 
itly) a single "effective viscosity" for the membrane [HIISHIES]- This approximation may not 
be as severe as it first appears. 




FIG. 3: Comparison of experimental (top) and simulated (bottom) phase separation kinetics for 
1:1:1 mixtures of DOPC/DPPC/Chol. Experimental figure adapted from Veatch and Keller (vesicle 
diameter is 30 /xm) See text for physical parameters and £ = 30 iJ,m, N = 1024 and At = 20/xs. 

The model uses eight physical parameters {T ,a,^,(t)Q,D^,rjf,rjm,C). The temperature T = 
21°C, box size £ = 30 fim and viscosity of water rjf = 0.01 Poise follow immediately from 
the experimental conditions. The remaining five parameters specify details of the specific 
bilayer system under study. Though precise measurements of all these parameters are not 
available, the numbers are known approximately, either for the specific DOPC/DPPC/Chol 
system under study or by analogy to other model systems expected to show similar behavior. 
00 = 0.4 may be determined from the DOPC/DPPC/Chol phase diagram [6]. The line 
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tension cr = 0.1 pN is based on flicker spectroscopy measurements of DOPC/DPPC/Chol 
mixtures [Ij, which also sets the rough order of the correlation length ,^ = 40 nm via the 
relation o"^ ^ ksT, a result motivated by the behavior of the Ising model [57] and believed 
to be consistent with experiments in ternary lipid systems |38]. The transport coefficients 
Vm = 5 X 10~^ surface Poise and = 7 x 10^^° cm^/s are difficult to precisely measure 
and are not known experimentally for the exact system studied here. The values chosen 
are consistent with measurements in related lipid systems [281 EH] and the requirement that 



The similarity between experiment and theory in Fig. |3] is striking and demonstrates 
that complex nonlinear bilayer dynamics over micron length scales and second time scales 
may be directly studied numerically. The displayed results depend upon both membrane and 
solvent viscosity; a naive 2D simulation assumes the wrong dissipation for length scales above 
Lsd (~ 1/im). More interestingly, dynamics of quasi-2D membrane phase separation are 
qualitatively different from the pure 2D case. We observe dynamical scaling with a continous 
morphology, with length scale R{t) ~ t^^^ for critical mixtures (1:1:1) at R{t) ^ Lgd, which 
can be explained using simple scaling arguments (as in [24j). In pure 2D hydrodynamic 
systems in the creeping flow limit, dynamical scale invariance breaks down [40j; in the 
absence of hydrodynamics, R{t) ~ t^^^ [211 [311 HQ]. Ri't) ~ ^^^^ scaling was also seen in [T5] : 
our results show that this scaling emerges naturally from quasi-2D hydrodynamics. Details 
of the unique scaling aspects of quasi-2D phase separation kinetics will be presented in a 
forthcoming paper. 

Our modeling is consistent with experimental [HI El EHl ESI EE] and theoretical results 
[HI HSl ESI ES] for a wide variety of phenomena in membrane biophysics. This approach lays 
the foundation for many extensions to more complex membrane and biomembrane systems. 
Our approach is readily combined with continuum simulations of out-of-plane membrane un- 
dulations [12] and coupling to the cytoskeleton pi] . The fluctuating- hydrodynamics scheme 
we have used is perfectly suited for an "immersed-boundary" [43j treatment of integral mem- 
brane proteins [44], in which stochastic fluid velocities (e.g. Eq. [2]) are directly coupled to 
protein motion. This sort of simulation may act as a bridge, using dynamics verifled in 
model membrane systems to elucidate our understanding of biomembranes. In particular, it 
will be interesting to further investigate both the thermodynamic nature [31] and dynamics 
of lipid raft models in light of the results obtained in this work. 
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Numerical details of simulation method 

Our overdamped model H for simulations of phase separation in a model membrane is 
given by 

SF 

(5^ + V . V)0(r, t) = MV'^j^ + e{r, t) (5) 



Vi{r,t) = j dPr' Tij{r - r') 



_5<p{v\t) 

where the continuum Fourier transform of Tjj is [HI [13] 



^1^V;.0(r',t) + O(r',t) 



(6) 



T,(,)^/.VT.(r)e-^^-r.-^^^ (7) 

where the integral is over all space. 

We follow Koga and Kawasaki [23] in using the fast Fourier transforms (FFT) as an 
efficient way to numerically evolve these dynamics. We modify their approach by using the 
Oseen tensor appropriate for the quasi-2D hydrodynamic environment as well as stochastic 
thermal forces. An C x C periodic geometry is assumed; the dynamics of the Fourier modes 
0q = d?r 0(r)e~*''''" follow from Eq. 6 

5t0q(t) + {v-V0(r,t)}q = -Afg2|^i^| (8) 

v^,{t) = T,,(q) |^^^V,0(r,t) + o| (9) 

{0^it)9*^'it')) = 2kBTMq^C\^^,Sit - t') (10) 
(Cq,(t)Q,,(t')) = 2kBTC^7]Uq' + q/Lsd)S^,6^,^'Sit - f). (11) 

where {/(r)}q is the Fourier transform of /(r) and * indicates complex conjugation. The 
variance of the Langevin forces 6 and ( are set by the fluctuation-dissipation theorem [211 12S| ■ 
These equations are solved numerically by truncating to N x N Fourier modes q = 
{'m,n)27c/C with —N/2 < m,n < N/2 (corresponding to a real space discretization size i = 
C/N). {. . .}q terms are evaluated in a hybrid real-space / Fourier-space fashion, handling 
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real-space derivatives and convolutions in q-space, local real-space operations in r-space and 
moving between the two representations via the FFT. 

Though we have written Eqs. [5]j6] as two separate equations, they only represent one 
dynamical equation, as the velocity field is set by the composition by Eq. [2j Substituting 
Eq. [6] into Eq. [5] yields a single Langevin equation for 0, but the coefficient of the thermal 
noise ( depends on through the V0 term of Eq. |5} This so-called "multiplicative noise" [26] 
should be treated via the Stratonovich interpretation, as it approximates a thermal force with 
a finite correlation time (set by the neglected fiuid inertia) [2S1 EH] • We use a semi- implicit 
Stratonovich integrator, which requires the nonconstant coefficient of the Langevin force to 
be averaged over its value at (p and an auxiliary value [30] ; the linear (but potentially most 
unstable) term is treated implicitly, and the nonlinear parts explicitly, as in semi-implicit 
solvers for the Cahn-Hilliard equation [31]. Our scheme is: 



Mt + At) = Mt) - ^ ^ M^^,^^ {^t {v'^* ■ V0(t) + v*^-- ■ (V0(t) + V0(t))/2 (12) 
-MV\r<p{t)-u<l){tf)}^-Q, 



<"W=7^..(q){^V,0} (13) 



1 



<r{t) = T.M)[-^/c,,j (14) 

kit) = ut) + - At {v*^-- ■ V0(t)}^ (15) 

(Bqe;) = 2kBTMq^C^At (16) 
= 2kBTC^ri^{q^ + g/L,,)At (17) 

Since 0(r) is a real field, we know that not all modes 0q are independent: 0* = </)_q. 
This means that (assuming is even), the modes (m,n) = (0, 0), (A^/2, 0), (0, A^/2), and 
{N/2, N/2) are guaranteed to be real. We choose these to be four of the required A^^ 
dynamical variables. The other independent modes are chosen to be (m,n) for — iV/2 < 
m < N/2 and < n < N/2, (m, 0) for < m < N/2, {m,N/2) for < m < N/2, and 
{N/2, n) for < n < N/2, as in [42J. The real and imaginary parts of each of these modes are 
both independent dynamic variables. The remaining modes are determined by the complex 
conjugates of the evolved modes. We also note that because the dynamics conserves total 
concentration, the mode (0, 0) must remain constant, and is not evolved. 

The random thermal forces 0{r,t) and (j{r,t) are also required to be real, which affects 
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their Fourier transforms, and therefore the variance of the integrals 0q(At) and Zqj(At). We 
know from the variance of 9q{t) (in the main paper) that (9q(At)6* (At)) = 2TMq^C'^At. 
If we write Bq = f + ig, we see that for the exphcitly real modes {m,n) = (0,0), (A^/2,0), 
(0, N/2), and (A^/2, N/2) where ^ = 0, = 2TMq'^C^At, but for complex modes, / and 

g are selected from a distribution with variance = (Is'P) = TMq^C'^At. The variance 

of is exactly analogous, but with (Zq^jZqj) = 2TC'^rjm{<f' + q/ Lsd)SijAt. 



The method of Eqs. p!2|fT7] allows stable time steps to be chosen that are nearly two 
orders of magnitude larger than those allowed by a simple explicit scheme, comparable to 
the gains seen in 
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